Finite-size corrections to the free energies of crystalline solids 
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We analyse the finite-size corrections to the free energy of crystals with a fixed center of mass. 
When we explicitly correct for the leading (lnA''/A*') corrections, the remaining free energy is found 
to depend linearly on 1/N. Extrapolating to the thermodynamic limit (A*" oo), we estimate 
the free energy of a defect free crystal of particles interacting through an r^^^-potential. We also 
estimate the free energy of perfect hard-sphere crystal near coexistence: at po-'^=1.0409, the excess 
free energy of a defect-free hard-sphere crystal is 5.91889(4)fcT per particle. This, however, is not 
the free energy of an equilibrium hard-sphere crystal. The presence of an finite concentration of 
vacancies results in a reduction of the free energy that is some two orders of magnitude larger than 
the present error estimate. 



The earliest numerical technique to compute the free- 
energy of crystalline solidsjiwas introduced some 30 years 
ago by Hoover and Reetlla. At present, the "single- 
occupancy-cell" method of Ree and Hoover is less widely 
used than the so-called "Einstein-crystal" method pro- 
posed by Frenkel and Ladcu. The latter method employs 
thermodynamic integration of the Helmholtz free energy 
along a reversible artificial pathway between the system 
of interest and an Einstein crystal. The Einstein crystals 
serves as a reference system, as its free energy can be 
computed analytically. Since its introduction, Einstein- 
crystal method has been used extensively in studies of 
phase equilibria involving crystalline solids. For numer- 
ical reasons - to suppress a weak divergence of the inte- 
grand - the Einstein-crystal method calculations have to 
be carried out at fixed center of mass. The free energy of 
the reference crystal is also calculated under the center- 
of-mass constraint, and the final calculated free energy 
of the unconstrained crystal is determined by correcting 
for the effect of imposing the constraint in the calcu- 
lations. In the original paper, the fixed center-of-mass 
constraint was only applied to the particle coordinates, 
but not to the corresponding momenta. This is irrele- 
vant as long as one computes the free-energy difference 
between two structures that have either both constrained 
or both unconstrained centers of mass. However, when 
computing the absolute free energy of a crystal, one needs 
to transform from the constrained to the unconstrained 
system. In the original paper, this transformation was 
not performed consistently. This resulted in a small but 
noticeable effect on the computed absolute free energy 
of the crystal. Below, we describe the proper approach 
to calculate the free energy of arbitrary molecular crys- 
talline solids. The derivation differs from the earlier work 



in two respects: first of all, we explicitly show the effect 
of momentum constraints. And secondly, we generalize 
the expression to an arbitrary crystal containing atoms 
or molecules with different masses. 

The main point of interest involves the calculation of 
the partition function of a crystal with and without a 
constrained center of mass. The partition function for 
an unconstrained, rf-dimensional crystalline solid of Nmoi 
molecules composed of a total of N atoms is given by 



Q^CN d^'^fd^'^peM-P^n^A)] 



(1) 



where CN^{h'"^"^'" Ni\N2l..N^\y^ , where there are Ni 
indistinguishable molecules of type 1, N2 molecules of 
type 2, etc., where A^i -I- + ■•• + Nm = Nmoi, and h 
is Planck's constant. It should be noted that, in all cal- 
culations of phase equilibria between systems that obey 
classical statistical mechanics, Planck's constant drops 
out of the result. Hence, in what follows, we omit all 
factors h. Using the result of the appendix in an article 
by Ryckaert and CiccottiQ, one can show that the con- 
strained partition function Q"^"" is given by 

Q™" = cw J d''''rd''''peM'mn,Pi)] 

X 5{a{f))5{G-^ ■ (2) 

where S{r) and a are the constraints and time deriva- 
tives of the constraints, respectively, and 



^ ruk dfk dfk 



k=l 



(3) 



The same integration limits implicit in Eqn. are also 
used in Eqn. (0). To constrain the center of mass (CM), 
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we take a{f) = J2iLi A**^*' ^"^^ thus, a = Y.i=iiP-i/'^i)&^ 
where Hi = rrii/ ^^rrii. Note that in Eqns. (|l|) and (g) 
we have assumed that there are no additional internal 
molecular constraints, such as fixed bond lengths or bond 
angles. 

We first consider the case of an Einstein crystal, 
which has a potential energy function given by U Em — 
(a/2) ^^-^(ri — ff^"')'^, where ff'"' are the equilibrium 
lattice positions. Note that the particles in a crystal are 
associated with specific lattice points and therefore be- 
have as if they are distinguishable - thus, cat = 1 (as we 
omit the factor l/h''-^^~^'>). It is easy to show that 



N 



nCM _ yCM pCM 



(4) 
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(5) 



and 



pCM 



N 



N 



i=l 1=1 



d/2 



Pi 



Ein ; 



(6) 



where M = ^ ^ , while ZEm and PEin are the corre- 
sponding contribution to QBimthe partition function of 
the unconstrained Einstein crystal. Clearly, 

Qgit! = {Y.m./J2^^r/\P'ahy4nr^'QEm (7) 



Similarly, one can show that the partition function for an 
arbitrary crystalline system subject to the CM constraint 
is given by 



N 



i=l 

with 

N 

Z""^' = J d'^^fexp[-/?f/(rl)]<5(^/i,rl) 



(8) 



(9) 



while the partition function of the unconstrained crystal 
is given by 



Q = Zl[{2nmJphY^^, 



with 



Z = J d'^'^fexp[-l3U{n)] 



(10) 



(11) 



Note that, as far as the kinetic part of the partition func- 
tion is concerned, the effect of the fixed center of mass 
constraint is the same for an Einstein crystal as for an 
arbitrary "realistic" crystal. Using Eqns. (^ and ([l^), 
the Helmholtz free energy difference between the con- 
strained and unconstrained crystal is given by 



/3(i^-F^*0 = -ln(0/Q''"'') 



CM\ 



We note that 

ZCM 



= ln(Z^*Vz) - - ln(27rM//3/i2) (12) 



/ dd^rexp[-(3U{ri)]d{J2^ ^^^n 



(13) 



where fcM = J^il^i^i^ ^^'^ 'P{'i'cm) is the probability 
distribution function of the center of mass, fcM- 

To calculate V[rcM) we exploit the fact that the equi- 
librium crystal lattice is invariant to translations over dis- 
placements through linear combinations of integer mul- 
tiples of the lattice vectors. This is true if the crystal 
lattice is subject to periodic boundary conditions. Conse- 
quently, the probability distribution of the center of mass 
of the lattice is evenly distributed over a volume equal 
to that of the Wigner-Seitz cell of the lattice positioned 
at the center of the volume over which we carry out the 
integration in the partition function. Since the average 
center of mass of the crystal is equal to the center of mass 
of the lattice, it follows that P{rcM) = 1/K;s = N^^s/V, 
where Vws is the volume of a Wigner-Seitz cell, and 
Nius is the number of such cells in the system. Thus, 
ZCM/Z = V{fcM = 0) = N.^s/V In the case of one 
molecule per cell, this implies Z'-^^'^ jZ — N,noi/V, where 
Nmoi is the number of molecules in the system. 

In the Frenkel-Ladd free energy calculation, the free 
energy difference between the constrained crystal and the 
reference systems is given by 



/3F' 



CM 



(14) 



where the statistical average of AU = UEin — U is calcu- 
lated by simulation for fixed CM as a function of A under 
an effective potential given by U{X) = (1 — X)U + XUEin- 
Note that the center of mass must be calculated in the 
same manner as described in the paragraph above. Fur- 
ther, note that this expression is only rigorously valid 
for systems interacting with continuous potentials. In 
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the case of particles with discontinuous potentials, e.g. 
hard particles, the internal potential energy cannot be 
turned off continuously. The calculation for this case 
differSpSlightly, and-.is discussed in detail in the original 
articled and in RcfB. 

Using Eqns. (0), ^ and (|l|), we find that the free 
energy per molecule of the unconstrained crystal is given 

by 



_0F_ 



|^)ln(2.//3«)-/^lnn 
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(3K^ 
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HV/N^oi) 



(15) 



If we consider the special case of a system of single- 
atom, identical particles (m^ — m and N — Nmoi), we 
obtain the following: 



1.97298(7)-5.17(6)/A' 
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/N + \n{N)/N vs. 1/N for a FCC 
r~^^) spheres at fcsT/e=1.0, and 



FIG. 1. f3F, 
crystal of soft (r ) spneres at kb 
per '^=1.1964, The solid line is a linear fit to the data. 
The coefficient of the 1/iV-term is: -5.17(6). 
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5.92 



d , , ^ dlniV 
ln(a/3/27r) 



where p = N/V. The difference between the present re- 
sult and the one obtained in ref.0 is in the fourth term 
on the right-hand side: —d\nN/2N. The original arti- 
cle implicitly gave the value +lnN/2N for a 3D crystal. 
While the difference between the two expressions tends 
to zero in the limit of large N, it is non-negligible for 
system sizes typically employed in the numerical calcu- 
lations. However, the calculated free energy differences 
between two solids, such as that between the FCC and 
HCP hard-sphere crystals, to which the method was ap- 
plied both in the original articleu and, more recently, in 
Rcf. ^, are unaffected by this correction. 

In practice, we usually need not calculate the abso- 
lute free energy of a crystal, but excess free energy, 
Fex = F — Fid, where Fid is the ideal gas free en- 
ergy. Let us therefore compute the finite-size correc- 
tions to the latter quantity: Given that (3Fid/N = 
-ln[y^(27rm//3/i2)''^/2/iV!]/Ar, we find 
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(17) 




where we have used In A^! sa TV In TV - + (In 2ttN) /2. 



FIG. 2. pFe^/N + \n{N)/N vs. 1/iV for an FCC 
crystal of hard spheres at a density /9cr^=l. 0409. The 
solid line is a linear fit to the data. The coefficient of 
the 1/Af-term is: -6.0(2). 

Hoover has analyzed the system-size dependence of 
the entropv. of a classical harmonic crystal with periodic 
boundariesQ. In this study, it was established that the 
leading finite-size correction to the free energy per parti- 
cle of a harmonic crystal is equal to /3~^ \nN/N. If the 
harmonic approximation is valid, then this implies that 
the integral in Eqn. ( p^ should vary as -|-lniV/7V plus 
higher order correction terms of the order of N~^, N~'^, 
etc. Consequently, an inspection of Eqn. ( p7| ) suggests 
that (3Fex/N+{d~-l) lnN/{2N) wiU scale as N~\ if we 
neglect terms of order 0{1/N'^). 

To test this prediction we have used the Einstcin- 
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crystal method to calculate the absolute Helmholtz free 
energy of a (three-dimensional) FCC crystal of soft 
spheres interacting with pair potential of u{r) = t{a /r)^'^ 
for systems of size 7V=216, 810, 1728, 5832 and 12096. 
The sizes were chosen such that the simulation box shape 
is cubic for each system. Further, the simulations were 
carried out at ksT / t=l.Q, /r?tT^ = 1.1964, and employed 
a coupling constant of a(T^/e=66.0. The results are 
shown in Figure ^ and are clearly consistent with the 
predictions. The solid line is a linear fit which extrap- 
olates to /3Fe:j,/A/'= 1.97298 (7) at = oo (where the 
figure between brackets is an estimate for the error in 
the last digit). Incidentally, we note that, at this den- 
sity and temperature, the FCC phase of soft-spheres, 
is more stable than the HCP phase (by an amount 
AFpcc-HCp/iNkBT) = 0.0028(8)). The present re- 
sults suggest that we are able to correctly account for 
the leading (In N/N) dependence of the free energy of an 
arbitrary crystal. In the analytical calculation of free- 
energy of a harmonic crystal, it is always assumed that 
the center of mass of the crystal is fixed. Hence, the 
numerical results presented above do not provide an in- 
dependent test of the validity of our expression for the 
contribution to the free energy due to the center-of-mass 
motion of the crystal. 

We can perform a similar analysis for a system of hard 
spheres (pcr^^ 1.0409^ X,nax = a/2=3000). The results 
are shown in Figure g. For hard spheres, f3Fex/N extrap- 
olates to a value of 5.91889(4) at iV = oo, well within the 
error margin of the original results of Hoover and ReeEl 
(5.924(15)). Note that the slopes of the fits (which are 
proportional to the 1/A^ behavior of the finite size ef- 
fect) are similar, although not exactly equal. It should 
be stressed that none of these calculations take into ac- 
count the existence of defects in the crystal, which, at 
these levels of precision, is significant. In fact, using 
the early numerical results by Bennett and AlderH we 
can estimate the equilibrium vacancy concentration in a 
hard-sphere crystal at coexistence to be 2.6 lO"**. Such 



a vacancy concentration has an noticeable effect on the 
location of the melting point. For instance, the Gibbs 
free energy per particle atucoexistence is lowered by an 
amount A/i sa —3 lO^^fcTa. This correction is far from 
negligible, as it is some two orders of magnitude larger 
than the present numerical accuracy in the absolute free 
energy. It is likely that vacancies also lower the equilib- 
rium free energy of the soft-sphere crystal. However, for 
that model, the equilibrium vacancy concentration has, 
to our knowledge, not been computed. 
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